clear all;

lambda_SOC=0.1*0.5;

spin_Delta=0.0;
spin_nvec=[0,0,1];
spin_nvec=spin_nvec/norm(spin_nvec);

sigx=[0,1;1,0];
sigy=[0,-i;i,0];
sigz=diag([1,-1]);
sig0=eye(2);

gen_SOC_term=@(bdvec) sqrt(2)*(i)*lambda_SOC*(bdvec(1)*sigx+bdvec(2)*sigy+bdvec(3)*sigz);


a1=[2.7508   -1.5882    4.4921];
a2=[0    3.1764    4.4921];
a3=[-2.7508   -1.5882    4.4921];

rec_a1=a1+a2-a3;
rec_a2=a1+a3-a2;
rec_a3=a3+a2-a1;

Rh_pos=[0.125000 0.125000 0.625000 ;
    0.625000 0.125000 0.125000 ;
    0.125000 0.625000 0.125000 ;
    0.125000 0.125000 0.125000];

Rh_pos=Rh_pos-0.125;

Rh_posC=Rh_pos;
Rh_posC(:,1)=Rh_pos(:,1)*a1(1)+Rh_pos(:,2)*a2(1)+Rh_pos(:,3)*a3(1);
Rh_posC(:,2)=Rh_pos(:,1)*a1(2)+Rh_pos(:,2)*a2(2)+Rh_pos(:,3)*a3(2);
Rh_posC(:,3)=Rh_pos(:,1)*a1(3)+Rh_pos(:,2)*a2(3)+Rh_pos(:,3)*a3(3);



Rh_posC



tetra1_posC=Rh_posC
tetra1_C=mean(tetra1_posC,1)

tetra2_posC=-Rh_posC
tetra2_C=mean(tetra2_posC,1)


newx_axis=zeros(4,3);
newy_axis=zeros(4,3);
newz_axis=zeros(4,3);
for ind=1:4
    zvec=tetra1_posC(ind,:)-tetra1_C;
    newz_axis(ind,:)=-zvec/norm(zvec);
end
newz_axis

newx_axis(4,:)=[1,0,0];
newx_axis(3,:)=[1,0,0];
newx_axis(2,:)=[-0.5,-0.5*sqrt(3),0];
newx_axis(1,:)=[-0.5,0.5*sqrt(3),0];


for ind=1:4
    xvec=newx_axis(ind,:);
    zvec=newz_axis(ind,:);
    newy_axis(ind,:)=cross(xvec,zvec);
end


dot(newy_axis(2,:),newz_axis(2,:))



Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(1,:))+[-1,0,1] %
Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(3,:))+[-1,1,0]
Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(4,:))+[-1,0,0]


Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(2,:))+[1,0,-1]
Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(3,:))+[0,1,-1]
Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(4,:))+[0,0,-1]

Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(1,:))+[0,-1,1]
Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(2,:))+[1,-1,0]
Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(4,:))+[0,-1,0]

Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(1,:))+[0,0,1]
Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(2,:))+[1,0,0]
Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(3,:))+[0,1,0]


tetra1=Rh_posC;
tetra1C=mean(tetra1,1);

tetra2=Rh_posC;
tetra2(1,:)=tetra2(1,:)-a3;
tetra2(2,:)=tetra2(2,:)-a1;
tetra2(3,:)=tetra2(3,:)-a2;

tetra1
tetra1C=mean(tetra1,1)
tetra2
tetra2C=mean(tetra2,1)




figure(1);
clf;


hold on;

linecc= [.5 .5 .5];
%plot_tetra(Rh_posC,0*a1,linecc)
% plot_tetra_v4(Rh_posC,a1,linecc,2)
% plot_tetra_v4(Rh_posC,a2,linecc,2)
% plot_tetra_v4(Rh_posC,a3,linecc,2)
% 
% %plot_tetra(-Rh_posC,0*a1,'r')
% %plot_tetra_v2(-Rh_posC,a1,linecc,2)
% %plot_tetra_v2(-Rh_posC,a2,linecc,2)
% %plot_tetra_v2(-Rh_posC,a3,linecc,2)
% 
% plot_tetra_v4(-Rh_posC,a1+a2,linecc,2)
% plot_tetra_v4(-Rh_posC,a2+a3,linecc,2)
% plot_tetra_v4(-Rh_posC,a1+a3,linecc,2)



hexpts1=[];
hexpts1(1,:)=Rh_posC(3,:)+a1;
hexpts1(2,:)=Rh_posC(2,:)+a2;
hexpts1(3,:)=Rh_posC(1,:)+a2;

hexpts1(4,:)=Rh_posC(3,:)+a3;
hexpts1(5,:)=Rh_posC(2,:)+a3;
hexpts1(6,:)=Rh_posC(1,:)+a1;


hexpts2=[];
hexpts2(1,:)=Rh_posC(4,:)+a1;
hexpts2(2,:)=Rh_posC(1,:)+a1;
hexpts2(3,:)=Rh_posC(2,:)+a3;
hexpts2(4,:)=Rh_posC(4,:)+a3;
hexpts2(5,:)=Rh_posC(1,:);
hexpts2(6,:)=Rh_posC(2,:);



hexpts3=[];
hexpts3(1,:)=Rh_posC(4,:)+a2;
hexpts3(2,:)=Rh_posC(2,:)+a2;
hexpts3(3,:)=Rh_posC(3,:)+a1;
hexpts3(4,:)=Rh_posC(4,:)+a1;
hexpts3(5,:)=Rh_posC(2,:);
hexpts3(6,:)=Rh_posC(3,:);



hexpts4=[];
hexpts4(1,:)=Rh_posC(4,:)+a3;
hexpts4(2,:)=Rh_posC(3,:)+a3;
hexpts4(3,:)=Rh_posC(1,:)+a2;
hexpts4(4,:)=Rh_posC(4,:)+a2;
hexpts4(5,:)=Rh_posC(3,:);
hexpts4(6,:)=Rh_posC(1,:);

hexpts1(7,:)=hexpts1(1,:);
for indh=1:6
    pt1=hexpts1(indh,:);
    pt2=hexpts1(indh+1,:);
    line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color','cyan','LineWidth',7); % [0,0,1]

end
%pp=fill3(hexpts1(:,1),hexpts1(:,2),hexpts1(:,3),[185,204,255]/255);pp.FaceAlpha=0.7;
pp=fill3(hexpts1(:,1),hexpts1(:,2),hexpts1(:,3),[0.4,0.4,1]);pp.FaceAlpha=0.5;


coorcenter=[-4.,-3.,hexpts1(1,3)-1];


coorcenter=[-2.,-5.,hexpts1(1,3)-1];
coordinatelen=1.5;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),coordinatelen,0,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,coordinatelen,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,0,coordinatelen,'MaxHeadSize',15,'Color','k','LineWidth',7)

axisvec = rec_a1/norm(rec_a1)*coordinatelen;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',15,'Color','k','LineWidth',7)
axisvec = rec_a2/norm(rec_a2)*coordinatelen;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',15,'Color','k','LineWidth',7)
axisvec = rec_a3/norm(rec_a3)*coordinatelen;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',15,'Color','k','LineWidth',7)




mirror_n1=(hexpts1(6,:)+hexpts1(1,:))/2-(hexpts1(3,:)+hexpts1(4,:))/2;
mirror_n1=mirror_n1/norm(mirror_n1);
mirror_n2=[0,0,1];

pltm1=4.;
pltm2=1.0;
pltmirror=zeros(4,3);
pltmirror(1,:)=-pltm1*mirror_n1-pltm2*mirror_n2;
pltmirror(2,:)=pltm1*mirror_n1-pltm2*mirror_n2;
pltmirror(3,:)=pltm1*mirror_n1+pltm2*mirror_n2;
pltmirror(4,:)=-pltm1*mirror_n1+pltm2*mirror_n2;
pltmirror(:,3)=pltmirror(:,3)+hexpts1(1,3);
%pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3;

mirror_line_color=[0,1,0,0.3];
mirror_line_color=[0.4660 0.6740 0.1880];
mirror_line_color=[0,0.5,0];

pt1=-pltm1*mirror_n1+[0,0,hexpts1(1,3)];
pt2=pltm1*mirror_n1+[0,0,hexpts1(1,3)];
pp=line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',mirror_line_color,'LineWidth',7)
%pp.FaceAlpha=0.3;

mirror_n1=(hexpts1(2,:)+hexpts1(1,:))/2-(hexpts1(5,:)+hexpts1(4,:))/2;
mirror_n1=mirror_n1/norm(mirror_n1);
mirror_n2=[0,0,1];

pltm1=4.;
pltm2=1.0;
pltmirror=zeros(4,3);
pltmirror(1,:)=-pltm1*mirror_n1-pltm2*mirror_n2;
pltmirror(2,:)=pltm1*mirror_n1-pltm2*mirror_n2;
pltmirror(3,:)=pltm1*mirror_n1+pltm2*mirror_n2;
pltmirror(4,:)=-pltm1*mirror_n1+pltm2*mirror_n2;
pltmirror(:,3)=pltmirror(:,3)+hexpts1(1,3);
%pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3;

pt1=-pltm1*mirror_n1+[0,0,hexpts1(1,3)];
pt2=pltm1*mirror_n1+[0,0,hexpts1(1,3)];
pp=line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',mirror_line_color,'LineWidth',7)
%pp.FaceAlpha=0.3;


mirror_n1=(hexpts1(2,:)+hexpts1(3,:))/2-(hexpts1(5,:)+hexpts1(6,:))/2;
mirror_n1=mirror_n1/norm(mirror_n1);
mirror_n2=[0,0,1];

pltm1=4.;
pltm2=1.0;
pltmirror=zeros(4,3);
pltmirror(1,:)=-pltm1*mirror_n1-pltm2*mirror_n2;
pltmirror(2,:)=pltm1*mirror_n1-pltm2*mirror_n2;
pltmirror(3,:)=pltm1*mirror_n1+pltm2*mirror_n2;
pltmirror(4,:)=-pltm1*mirror_n1+pltm2*mirror_n2;
pltmirror(:,3)=pltmirror(:,3)+hexpts1(1,3);
%pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3;

pt1=-pltm1*mirror_n1+[0,0,hexpts1(1,3)];
pt2=pltm1*mirror_n1+[0,0,hexpts1(1,3)];
pp=line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',mirror_line_color,'LineWidth',7)
%pp.FaceAlpha=0.3;






sig_x=[0,1;1,0];
sig_y=[0,-i;i,0];
sig_z=[1,0;0,-1];
sig_0=eye(2);

gen_spinh=@(nnvec) (sig_x*nnvec(1)+sig_y*nnvec(2)+sig_z*nnvec(3))/norm(nnvec);

eva_spinS =@(state_wf) [state_wf'*sig_x*state_wf,state_wf'*sig_y*state_wf,state_wf'*sig_z*state_wf];



newnvec_tetra=Rh_posC(3,:)-tetra1C;
newnvec_tetra=newnvec_tetra/norm(newnvec_tetra);





tripts=Rh_posC(1:3,1:3);
tripts_xaxis=newx_axis(1:3,1:3);
%quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=1:3
    %plt_porbital(tripts(ind,:),tripts_xaxis(ind,:));
end





shiftvec=a2;
tripts=Rh_posC([1,2,4],1:3);
%tripts_xaxis=newx_axis([1,2,4],1:3);

newnvec=Rh_posC(3,:)-tetra1_C;
newnvec=newnvec/norm(newnvec);
newnvecx=[1,0,0];
newnvecy=cross(newnvec,newnvecx);
tripts_xaxis=zeros(3,3);
tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
tripts_xaxis(3,:)=[-1,0,0];

tripts_zaxis=zeros(3,3);
for ind=1:3
    tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;

end

for ind=[1,2]
    %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
    plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);

end

%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')


R3mat=zeros(3);
R3mat(3,3)=1;
R3mat(1,1)=-0.5;
R3mat(2,2)=-0.5;
R3mat(1,2)=0.5*sqrt(3);
R3mat(2,1)=-0.5*sqrt(3);


tripts_xaxis=(R3mat*tripts_xaxis')';
shiftvec=a1;
tripts=Rh_posC([3,1,4],1:3);
%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')

tripts_zaxis=zeros(3,3);
for ind=1:3
    tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;

end



for ind=1:3
    %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
end

for ind=[1,2]
    %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
    plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);

end


tripts_xaxis=(R3mat*tripts_xaxis')';
shiftvec=a3;
tripts=Rh_posC([2,3,4],1:3);
%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=1:3
    %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
end


tripts_zaxis=zeros(3,3);
for ind=1:3
    tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;

end

for ind=[1,2]
    %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
    plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);

end


%plt_D2_orbital([0,0,0],[0,0,0]);

'add light'

%set(gca, 'FontSize', 10, 'XColor', 'black', 'YColor', 'black', 'ZColor', 'black', 'LineWidth', 1);
%axis equal
%lightangle(-45,30);
%campos([-1.0926  -1.7493  0.9397]);

%  camlight('headlight')







% plt_D0_orbital([0,0,1],[1,0,0],[0,1,1],-1);

%view([1,1,0.5])

%view([139,22])
%view([0,0,1])

view([20,35])

set(gca,'visible','off')


material dull
 camlight('headlight')

hold off;


box on;
axis equal;




[az,el] = view



%%


set(gca,'visible','off')
set(gcf, 'color', 'none');
set(gca, 'color', 'none');


exportgraphics(gcf,'Fig4c-pzCLS-hex.eps','ContentType','image','BackgroundColor','none','Resolution',400)
% 'ContentType','vector'


%%

figure(2);

clf;

hold on;


coorcenter=[-2.,-5.,hexpts1(1,3)-1];
coordinatelen=1.5;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),coordinatelen,0,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,coordinatelen,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,0,coordinatelen,'MaxHeadSize',15,'Color','k','LineWidth',7)

axisvec = rec_a1/norm(rec_a1)*coordinatelen;
quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',40,'Color','k','LineWidth',15)
axisvec = rec_a2/norm(rec_a2)*coordinatelen;
quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',40,'Color','k','LineWidth',15)
axisvec = rec_a3/norm(rec_a3)*coordinatelen;
quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',40,'Color','k','LineWidth',15)

view([20,35])

set(gca,'visible','off')


material dull
 camlight('headlight')

hold off;


box on;
axis equal;



%%


set(gca,'visible','off')
set(gcf, 'color', 'none');
set(gca, 'color', 'none');


exportgraphics(gcf,'Fig4c-pzCLS-hex-coordinate-frame.eps','ContentType','vector','BackgroundColor','none','Resolution',400)
% 'ContentType','vector'




%%

% %%
% 
% 
% %exportgraphics(gcf,'Pyrochlore_dz2_loopstate.eps','ContentType','vector','BackgroundColor','none')
% 
% figure(1)
% clf
% 
% hold on;
% 
% 
% linecc= [.5 .5 .5];
% %plot_tetra(Rh_posC,0*a1,linecc)
% plot_tetra_v3(Rh_posC,a2,linecc,2)
% 
% 
% shiftvec=a2;
% tripts=Rh_posC([1,2,4],1:3);
% %tripts_xaxis=newx_axis([1,2,4],1:3);
% 
% newnvec=Rh_posC(3,:)-tetra1_C;
% newnvec=newnvec/norm(newnvec);
% newnvecx=[1,0,0];
% newnvecy=cross(newnvec,newnvecx);
% tripts_xaxis=zeros(3,3);
% tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
% tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
% tripts_xaxis(3,:)=[-1,0,0];
% 
% tripts_zaxis=zeros(3,3);
% for ind=1:3
%     tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;
% 
% end
% 
% for ind=[1,2]
%     %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
%     plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);
% 
% end
% ind=3
% plt_D0_orbital_v2(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);
% 
% 
% pltmirror=[0,-1,-1.2;
%     0,-1,0.6;
%     0,1,0.6;
%     0,1,-1.2]*2.5;
% 
% pltmirror(:,1)=pltmirror(:,1)+tetra1_C(1)+shiftvec(1);
% pltmirror(:,2)=pltmirror(:,2)+tetra1_C(2)+shiftvec(2);
% pltmirror(:,3)=pltmirror(:,3)+tetra1_C(3)+shiftvec(3);
% 
% 
% pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3;
% 
% 
% 
% 
% pt1=tripts(1,:)+shiftvec;
% pt2=tripts(2,:)+shiftvec;
% line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color','k','LineWidth',7)
% 
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% 
% 
% R3mat=zeros(3);
% R3mat(3,3)=1;
% R3mat(1,1)=-0.5;
% R3mat(2,2)=-0.5;
% R3mat(1,2)=0.5*sqrt(3);
% R3mat(2,1)=-0.5*sqrt(3);
% 
% 
% 
% 
% view([0,0,1])
% 
% set(gca,'visible','off')
% 
% %camlight('headlight')
% 
% 
% 
% hold off;
% 
% 
% box on;
% axis equal;
% 
% [az,el] = view
